Example Single Cell RNA-seq Analysis

Author

Kazi Tanvir Hasan

Published

March 19, 2023

Attaching SeuratObject

This is an example analysis using Seurat package to analyze scRNA-seq data in R. Seurat is a popular package for single-cell RNA sequencing data analysis that provides a wide range of functions for preprocessing, normalization, quality control, clustering, and visualization of scRNA-seq data.

Import Seurat Object

# Load the Seurat dbject dataset
objectSeurat <- readRDS(
  "./data/AD00201.rds"
)

Data QC and Inspection

To perform quality control on the data after importing (or creating) the Seurat object, follow these steps:

  1. Check for cells with too few genes detected. These cells may not have been sequenced deeply enough for reliable characterization. Filter them out.

  2. Check for cells with too many genes detected. These cells may represent doublets or multiplets, where two or more cells are in the same droplet, and therefore share the same cell barcode. Filter them out.

  3. Check for cells with a high mitochondrial transcript percentage. As scRNA-seq experiments use oligo-T to capture mRNAs, mitochondrial transcripts should be relatively under-representative due to their lack of poly-A tails. However, some mitochondrial transcripts may be captured. Cells with high mitochondrial transcript percentages likely represent cells under stress (e.g., hypoxia) that produce a lot of mitochondria or an abnormally high amount of truncated mitochondrial transcripts. Filter them out.

Note that Seurat automatically summarizes the numbers of detected genes when creating the Seurat object (nFeature_RNA) and the number of detected transcripts (nCount_RNA). However, you need to calculate mitochondrial transcript percentages manually. Seurat provides an easy solution.

# The [[ operator can add columns to object metadata. 
# This is a great place to stash QC stats
objectSeurat[["percent.mt"]] <- PercentageFeatureSet(
  objectSeurat, pattern = "^MT-"
)

When filtering cells based on quality control (QC) metrics, it is important to note that there is no one-size-fits-all filtering criteria. The normal ranges of these metrics can vary dramatically from one experiment to another, depending on sample origin, reagents, and sequencing depths. One suggestion to address this variability is to only filter out outlier cells, which are the minority of cells with certain QC metrics clearly above or below the majority of cells.

To identify these outlier cells, we first need to know how the QC metric values are distributed in the data. We can visualize this distribution using a violin plot for each of the metrics. The violin plot displays the distribution of values for each metric, with the width of the plot showing the density of cells at each value. By creating these plots, we can identify cells that fall far outside the main distribution of values for each metric, indicating that they are likely outliers. These outlier cells can then be filtered out, while retaining the remaining cells for downstream analyses.

# Visualize QC metrics as a violin plot
VlnPlot(
  objectSeurat, 
  features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
  ncol = 3
)

Or if you don’t like the dots (individual cells)

VlnPlot(
  objectSeurat, 
  features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
  ncol = 3, 
  pt.size = 0
)

Based on our prior knowledge of these metrics, we would expect the number of detected genes and number of detected transcripts to be well correlated, while mitochondrial transcript percentage is not expected to be strongly correlated with the other two metrics.

# FeatureScatter is typically used to visualize feature-feature relationships,
# but can be used for anything calculated by the object, i.e. columns in object
# metadata, PC scores etc.

plot1 <- FeatureScatter(
  objectSeurat, 
  feature1 = "nCount_RNA", 
  feature2 = "percent.mt"
)

plot2 <- FeatureScatter(
  objectSeurat, 
  feature1 = "nCount_RNA", 
  feature2 = "nFeature_RNA"
)

plot1 + plot2

The R package patchwork was created to simplify the organization of plots generated by ggplot2, which Seurat relies on for its own plotting functions. Without patchwork, attempting to execute plot1 + plot2 would be invalid.

To perform quality control (QC) on our data, we only need to establish a cutoff for one of two metrics: gene count or transcript count. Additionally, we should set an upper threshold for the percentage of mitochondrial transcripts. For this dataset, appropriate parameters could be a gene count between 500 and 5000, and a mitochondrial transcript percentage below 5%, though other thresholds may also be effective.

objectSeurat <- subset(
  objectSeurat, 
  subset = nFeature_RNA > 500 & nFeature_RNA < 5000 & percent.mt < 5
)

It is important to note that in some cases, additional QC measures may be necessary. One possible concern is the occurrence of doublets, which may not exhibit a higher count of detected genes or transcripts due to variations in captured RNA across cells. To address this, several tools have been developed to identify potential doublets or multiplets. For example, Doublet Finder creates artificial doublets by averaging cells at random and compares each cell to these artificial doublets to determine whether it is more similar to them or not. This approach can help determine whether a cell is likely a doublet. Furthermore, filtering based on mitochondrial transcript percentage alone may not suffice to eliminate unhealthy or stressed cells. In such cases, additional filtering techniques, such as those utilizing machine learning-based predictions, may be necessary.

Normalization and Data scaling

Normalization

Similar to bulk RNA-seq, the amount of RNA captured can vary from cell to cell, which means that comparing the number of captured transcripts for each gene directly between cells is not appropriate. To enable the comparison of gene expression levels across different cells, a normalization step is necessary. The most commonly used normalization method in scRNA-seq data analysis is quite similar to TPM (Transcripts Per Million reads). It involves normalizing the feature expression measurements for each cell to the total expression and then multiplying this by a scale factor (typically 10000). Finally, the resulting expression levels are log-transformed to better fit a normal distribution. Notably, before log-transforming the data, a pseudocount is added to each value to ensure that genes with zero transcripts detected in a cell still have values of zero after log-transformation.

objectSeurat <- NormalizeData(
  objectSeurat, 
  normalization.method = "LogNormalize"
)

Feature selection

Single-cell RNA sequencing (scRNA-seq) has a significant advantage over bulk RNA-seq as it allows for the identification of cellular heterogeneity by looking for cell groups with distinct molecular signatures. However, not all genes are equally informative or contribute equally to the identification of different cell groups. Genes with low expression levels or those with similar expression levels across all cells may not provide much information and can dilute differences between distinct cell groups. Therefore, before exploring scRNA-seq data further, it is essential to perform a proper feature selection process.

In Seurat and other scRNA-seq data analysis methods, this process involves the identification of highly variable features or genes. These are genes with the most varied expression levels across cells and are used to distinguish different cell groups with distinct molecular signatures.

objectSeurat <- FindVariableFeatures(
  objectSeurat, 
  selection.method = "vst", 
  nfeatures = 2000
)
#https://rdrr.io/bioc/DESeq2/man/varianceStabilizingTransformation.html

The VST function is used to transform the count data in scRNA-seq data analysis. This transformation aims to make the data approximately homoskedastic, meaning that the variance is constant across the range of mean values. It also normalizes the data with respect to library size. However, the transformation can be sensitive to size factors, which can vary widely. To address this issue, Seurat provides another transformation called rlog, which is less sensitive to size factors. Both transformations can be useful when checking for outliers or when using machine learning techniques such as clustering or linear discriminant analysis.

In Seurat, the VST function is used to identify highly variable features/genes by calculating the standardized variance of each gene across cells and picking the top ones as highly variable features. By default, the top 2000 highly variable features are selected, but this number can be changed using the nfeatures option.

There is no good criteria to determine how many highly variable features to use. The number of highly variable features to use is determined through iteration and observation of the results. Typically, a value between 2000 to 5000 is sufficient and using a different value does not significantly affect the results.

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(objectSeurat), 10)

top10
 [1] "CTGF"         "CTC-340A15.2" "NPY"          "KIF26B"       "SRGN"        
 [6] "CCL3"         "CNR1"         "RELN"         "AJ006998.2"   "AC011288.2"  
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(objectSeurat)

plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
When using repel, set xnudge and ynudge to 0 for optimal results

plot1
Warning: Transformation introduced infinite values in continuous x-axis


plot2
Warning: Transformation introduced infinite values in continuous x-axis

Data Scaling

When analyzing single-cell RNA-seq data, it is important to account for differences in gene expression levels and distributions. Without doing so, highly expressed genes would dominate the analysis, potentially obscuring important information from lower expressed genes. To prevent this, a scaling transformation is applied to the data using a set of selected features. This scaling allows for a fairer comparison between genes, ensuring that each gene’s contribution to the analysis is given appropriate weight. Scaling is a common technique used across many different data science fields to ensure that all features are given equal consideration.

all.genes <- rownames(objectSeurat)
objectSeurat <- ScaleData(objectSeurat, features = all.genes)
Centering and scaling data matrix

Data Clustering (PCA/UMAP)

After identifying highly variable genes and scaling the scRNA-seq data, it is strongly recommended to apply linear dimension reduction before conducting further analysis. This step has several benefits, including but not limited to:

  1. Compressing the data and making it much more compact, which results in faster computation times.

  2. Enhancing signal robustness by summarizing measurements of related features, which is especially important as scRNA-seq data is intrinsically sparse.

objectSeurat <- RunPCA(
  objectSeurat, 
  features = VariableFeatures(object = objectSeurat)
)
PC_ 1 
Positive:  SYT1, RASGEF1B, GRIN1, RBFOX1, MT-CO1, MEG3, MT-ND1, MT-CO3, SLC26A3, DLGAP1 
       SNAP25, STXBP5L, CACNA1B, BCYRN1, SYN2, INO80D, PHYHIP, MT-CO2, CSMD1, DNM1 
       LINC01609, ENC1, MT-ATP6, SLC17A7, KCNQ5, WBSCR17, IQCJ-SCHIP1, CHD5, MT-ND2, RYR2 
Negative:  PLP1, SLC1A3, ST18, APOE, CST3, ERBB4, NEAT1, MIR219A2, MT2A, APOD 
       SLC1A2, PREX2, RNF219-AS1, B2M, ADGRV1, LINC01608, MOBP, SPP1, ATP1A2, LSAMP-AS1 
       CNDP1, FGFR3, AQP4, NHSL1, PTPRZ1, LINC00499, EMX2, COL5A3, ETNPPL, AC002429.5 
PC_ 2 
Positive:  PTPRD, ST18, APOD, PLP1, LHFPL3, CNTNAP2, MIR219A2, LRRTM4, NXPH1, OPCML 
       PCDH15, MOBP, LINC01608, GRIK1, SGCZ, LUZP2, RP4-668E10.4, MMP16, ADARB2, SEMA5A 
       GRIA4, PDE1A, CA10, LINC01170, SMOC1, CNDP1, GRIK2, SPP1, KAZN, TNR 
Negative:  ADGRV1, RNF219-AS1, GPC5, SLC1A2, TPD52L1, LINC00499, AQP4, SLC4A4, BMPR1B, EMX2 
       COL5A3, SLC1A3, TRPM3, NKAIN3, GLIS3, ETNPPL, PAMR1, FGFR3, NHSL1, MT2A 
       STON2, APOE, CLU, AC002429.5, RANBP3L, GNA14, MGST1, GJA1, ATP1A2, SLCO1C1 
PC_ 3 
Positive:  PTPRZ1, NXPH1, PCDH15, SOX6, LUZP2, GRIK1, LHFPL3, VCAN, SGCZ, MMP16 
       RP4-668E10.4, BRINP3, ERBB4, GRIK2, OPCML, TNR, LINC01322, KCNMB2-AS1, NKAIN3, SEMA5A 
       KAZN, XKR4, KIAA1211, GPC6, SMOC1, CA10, COL11A1, CSMD1, ADARB2, MIR3681HG 
Negative:  ST18, PLP1, MOBP, NKAIN2, LINC01608, SPP1, MIR219A2, LPAR6, ADAM28, PDE1A 
       CNDP1, DOCK8, FYB, APBB1IP, NPTX2, CD74, BDNF, MAML3, SRGN, NPAS4 
       PTPRC, EGR4, INPP5D, SYK, BLNK, HLA-DRA, OLR1, CX3CR1, ATP8B4, SAMSN1 
PC_ 4 
Positive:  GAD2, ARL4C, GAD1, LINC00486, IGF1, ADARB2, GRIK1, FOS, RP11-123O10.4, DOCK8 
       ERBB4, DLX6-AS1, BTBD11, KCNIP1, MTRNR2L8, ADAM28, ZNF385D, GRIP1, DLX1, SDK1 
       SCG2, EGR1, LPAR6, APBB1IP, CD74, SOX1, FOSB, FYB, RP11-154D17.1, LINC01609 
Negative:  AC011288.2, LDB2, PTPRD, RALYL, ST6GALNAC5, NKAIN2, IQCJ-SCHIP1, RP11-380P13.1, LY86-AS1, FAM19A1 
       SV2B, AC114765.1, PDE1A, LINC01250, MLIP, KCNIP4, ANO3, CBLN2, CADPS2, OLFM3 
       FSTL4, SLIT3, CHN1, MKL2, RP11-191L9.4, AC067956.1, AJ006998.2, MIR137HG, HECW1, PHACTR1 
PC_ 5 
Positive:  DOCK8, ADAM28, ST6GAL1, LPAR6, FYB, APBB1IP, CD74, SRGN, SYK, SRGAP2B 
       C10orf11, INPP5D, PTPRC, RP11-624C23.1, SRGAP2, CSF2RA, HLA-DRA, HS3ST4, P2RY12, BLNK 
       OLR1, ATP8B4, ARHGAP15, SAMSN1, A2M, CX3CR1, TBXAS1, CSF3R, CD86, RP11-480C22.1 
Negative:  PLP1, ST18, MIR219A2, MOBP, ARC, LINC01608, FOSB, NPAS4, EGR4, MTRNR2L8 
       CNDP1, EGR1, ERBB4, JUND, NPTX2, EGR3, PIM3, KIRREL3, LINC01170, DUSP5 
       INHBA, NKAIN2, RP11-154D17.1, APOD, RP11-81H3.2, SLC27A6, MIDN, KDM6B, BCOR, LBH 

# Examine and visualize PCA results a few different ways
print(objectSeurat[["pca"]], dims = 1:2, nfeatures = 5)
PC_ 1 
Positive:  SYT1, RASGEF1B, GRIN1, RBFOX1, MT-CO1 
Negative:  PLP1, SLC1A3, ST18, APOE, CST3 
PC_ 2 
Positive:  PTPRD, ST18, APOD, PLP1, LHFPL3 
Negative:  ADGRV1, RNF219-AS1, GPC5, SLC1A2, TPD52L1 
VizDimLoadings(
  objectSeurat, dims = 1:2, 
  nfeatures = 15, 
  reduction = "pca"
)

DimPlot(objectSeurat, reduction = "pca")


DimHeatmap(objectSeurat, dims = 1, cells = 500, balanced = TRUE)

When analyzing scRNA-seq data, it is important to identify the underlying systematic patterns of variation, which can be captured by latent variables such as principal component analysis (PCA) or factor analysis (FA). However, it’s not enough to simply identify these patterns - we also need to test for associations between observed variables and these latent variables. This is where the jackstraw method comes in - it enables us to statistically test for such associations, using PCs or other estimates as our latent variables.

# NOTE: This process can take a long time for big datasets, comment out for expediency. More
# approximate techniques such as those implemented in ElbowPlot() can be used to reduce
# computation time
objectSeurat <- JackStraw(objectSeurat, num.replicate = 100)

objectSeurat <- ScoreJackStraw(objectSeurat, dims = 1:20)

JackStrawPlot(objectSeurat, dims = 1:15)
Warning: Removed 21000 rows containing missing values (`geom_point()`).

ElbowPlot(objectSeurat)

In our analysis, we selected 10 principal components (PCs) to capture the systematic variation in the data. However, we recommend users to test downstream analyses using different numbers of PCs, such as 10, 15, or 50, as we have observed that the results do not differ significantly. It is important to note that choosing a lower number of PCs, such as 5, can negatively impact the downstream analysis results. Hence, we suggest erring on the higher side when choosing the number of PCs to be used.

objectSeurat <- FindNeighbors(objectSeurat, dims = 1:10)
Computing nearest neighbor graph
Computing SNN
objectSeurat <- FindClusters(objectSeurat, resolution = 0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 8970
Number of edges: 296870

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9382
Number of communities: 15
Elapsed time: 0 seconds
# Look at cluster IDs of the first 5 cells
head(Idents(objectSeurat), 5)
AAACCTGAGTCACGCC_Ct5_Ct6 AAACCTGCATGGATGG_Ct5_Ct6 AAACCTGGTAGCTAAA_Ct5_Ct6 
                       1                        0                        0 
AAACCTGGTGCTAGCC_Ct5_Ct6 AAACCTGGTGTCGCTG_Ct5_Ct6 
                      13                        1 
Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14

When performing linear dimension reduction, there are both advantages and disadvantages to consider. One of the benefits is that each Principal Component (PC) is a linear combination of gene expressions, making the interpretation of PCs straightforward. Additionally, the data is compressed without distorting it, thus preserving most of the information in the data. However, a potential drawback is that more than 10 PCs are typically needed to capture most of the information. While this is generally acceptable for most analyses, it can be problematic for visualization since ordinary people can only perceive up to three dimensions.

Non-linear dimension reduction techniques such as t-distributed Stochastic Neighbor Embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) are commonly used in scRNA-seq data analysis to address this issue. In this example analysis, only the umap method is used.

objectSeurat <- RunUMAP(objectSeurat, dims = 1:10)
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
19:58:37 UMAP embedding parameters a = 0.9922 b = 1.112
19:58:37 Read 8970 rows and found 10 numeric columns
19:58:37 Using Annoy for neighbor search, n_neighbors = 30
19:58:37 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
19:58:37 Writing NN index file to temp file /var/folders/xn/zm5z4n5j6t5fbn1m6pzwtqsc0000gn/T//RtmpPtXEx8/file12b1713ef224
19:58:37 Searching Annoy index using 1 thread, search_k = 3000
19:58:39 Annoy recall = 100%
19:58:40 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
19:58:40 Initializing from normalized Laplacian + noise (using irlba)
19:58:40 Commencing optimization for 500 epochs, with 369150 positive edges
19:58:53 Optimization finished

DimPlot(objectSeurat, reduction = "umap")

Markers Identification

# find all markers of cluster 1
cluster1.markers <- FindMarkers(
  objectSeurat, 
  ident.1 = 1, 
  min.pct = 0.25
)
For a more efficient implementation of the Wilcoxon Rank Sum Test,
(default method for FindMarkers) please install the limma package
--------------------------------------------
install.packages('BiocManager')
BiocManager::install('limma')
--------------------------------------------
After installation of limma, Seurat will automatically use the more 
efficient implementation (no further action necessary).
This message will be shown once per session

head(cluster1.markers, n = 5)
        p_val avg_log2FC pct.1 pct.2 p_val_adj
PLPP3       0       1.65 0.512 0.098         0
NFIA        0       1.50 0.874 0.360         0
CACHD1      0       1.92 0.556 0.106         0
C1orf61     0       2.00 0.805 0.227         0
ATP1A2      0       2.27 0.686 0.086         0

VlnPlot(
  objectSeurat, 
  features = c(
    row.names(cluster1.markers)[1],
    row.names(cluster1.markers)[2]
  )
)

# find all markers of cluster 2
cluster2.markers <- FindMarkers(
  objectSeurat, 
  ident.1 = 2,
  min.pct = 0.25
)

head(cluster2.markers, n = 5)
        p_val avg_log2FC pct.1 pct.2 p_val_adj
STXBP5L     0       1.96 0.886 0.371         0
KCNIP4      0       1.87 0.969 0.511         0
KCNQ5       0       1.99 0.761 0.263         0
SYT1        0       1.79 0.957 0.398         0
MEG3        0       1.67 0.992 0.605         0

VlnPlot(
  objectSeurat, 
  features = c(
    row.names(cluster2.markers)[1], 
    row.names(cluster2.markers)[2]
  )
)

# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(
  objectSeurat, 
  ident.1 = 5, ident.2 = c(0, 3), 
  min.pct = 0.25
)

head(cluster5.markers, n = 5)
         p_val avg_log2FC pct.1 pct.2 p_val_adj
CHRM3        0       3.02 0.765 0.051         0
INO80D       0       3.83 0.663 0.058         0
SYNPR        0       3.60 0.765 0.028         0
ROBO2        0       4.23 0.925 0.093         0
RASGEF1B     0       5.71 0.683 0.036         0

VlnPlot(
  objectSeurat, 
  features = c(
    row.names(cluster5.markers)[1], 
    row.names(cluster5.markers)[2]
  )
)

# find markers for every cluster compared to all remaining cells, report only the positive ones
objectSeurat.markers <- FindAllMarkers(
  objectSeurat,
  only.pos = TRUE, 
  min.pct = 0.25, 
  logfc.threshold = 0.25
)
Calculating cluster 0
Calculating cluster 1
Calculating cluster 2
Calculating cluster 3
Calculating cluster 4
Calculating cluster 5
Calculating cluster 6
Calculating cluster 7
Calculating cluster 8
Calculating cluster 9
Calculating cluster 10
Calculating cluster 11
Calculating cluster 12
Calculating cluster 13
Calculating cluster 14
x <- objectSeurat.markers %>% 
  group_by(cluster) %>%
  top_n(n = 1, wt = avg_log2FC)

FeaturePlot(objectSeurat, features = x$gene[1:4])


FeaturePlot(objectSeurat, features = x$gene[5:8])

p <- FeaturePlot(
  objectSeurat,
  features = c(
    "MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ",
    "PPBP", "CD8A"
  ),
  combine = FALSE
)
Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
: The following requested variables were not found: MS4A1, FCER1A, PPBP

p <- lapply(X = p, FUN = function(x) x +
              theme(plot.title = element_text(size = 8)) +
              theme(axis.title.y = element_text(size = 5)) +
              theme(axis.title.x = element_text(size = 5)) +
              theme(axis.text.y = element_text(size = 5)) +
              theme(axis.text.x = element_text(size = 5)) +
              theme(legend.position = "none")  )

CombinePlots(plots = p)
Warning: CombinePlots is being deprecated. Plots should now be combined using
the patchwork system.

top10 <- objectSeurat.markers %>% 
            group_by(cluster) %>% 
            top_n(n = 10, wt = avg_log2FC)

top10
# A tibble: 150 × 7
# Groups:   cluster [15]
   p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene     
   <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>    
 1     0       4.31 0.991 0.219         0 0       PLP1     
 2     0       3.98 0.911 0.107         0 0       ST18     
 3     0       3.72 0.963 0.285         0 0       CTNNA3   
 4     0       3.50 0.962 0.271         0 0       MBP      
 5     0       3.49 0.853 0.108         0 0       TMEM144  
 6     0       3.25 0.609 0.045         0 0       LINC01608
 7     0       3.19 0.748 0.113         0 0       MIR219A2 
 8     0       3.13 0.633 0.059         0 0       MOBP     
 9     0       3.09 0.681 0.092         0 0       TF       
10     0       3.08 0.696 0.1           0 0       CLDN11   
# … with 140 more rows
p2 <- DoHeatmap(
  objectSeurat,
  features = top10$gene, 
  group.bar.height = 0.01,
  size = 3,
  combine = FALSE
) 

p2 <- lapply(X = p2, FUN = function(x) x + 
               theme(plot.title = element_text(size = 8)) +
               theme(axis.title.y = element_text(size = 5)) +
               theme(axis.title.x = element_text(size = 5)) +
               theme(axis.text.y = element_text(size = 3)) +
               theme(legend.position = "none")  )

CombinePlots(plots = p2)
Warning: CombinePlots is being deprecated. Plots should now be combined using
the patchwork system.

sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lubridate_1.9.2    forcats_1.0.0      stringr_1.5.0      dplyr_1.1.0       
 [5] purrr_1.0.1        readr_2.1.4        tidyr_1.3.0        tibble_3.2.0      
 [9] ggplot2_3.4.1      tidyverse_2.0.0    patchwork_1.1.2    SeuratObject_4.1.3
[13] Seurat_4.3.0       conflicted_1.2.0  

loaded via a namespace (and not attached):
  [1] Rtsne_0.16             colorspace_2.1-0       deldir_1.0-6          
  [4] ellipsis_0.3.2         ggridges_0.5.4         rstudioapi_0.14       
  [7] spatstat.data_3.0-1    farver_2.1.1           leiden_0.4.3          
 [10] listenv_0.9.0          ggrepel_0.9.3          fansi_1.0.4           
 [13] codetools_0.2-19       splines_4.2.1          cachem_1.0.7          
 [16] knitr_1.42             polyclip_1.10-4        jsonlite_1.8.4        
 [19] ica_1.0-3              cluster_2.1.4          png_0.1-8             
 [22] uwot_0.1.14            shiny_1.7.4            sctransform_0.3.5     
 [25] spatstat.sparse_3.0-1  compiler_4.2.1         httr_1.4.5            
 [28] Matrix_1.5-3           fastmap_1.1.1          lazyeval_0.2.2        
 [31] cli_3.6.0              later_1.3.0            htmltools_0.5.4       
 [34] tools_4.2.1            igraph_1.4.1           gtable_0.3.2          
 [37] glue_1.6.2             RANN_2.6.1             reshape2_1.4.4        
 [40] Rcpp_1.0.10            scattermore_0.8        vctrs_0.5.2           
 [43] nlme_3.1-162           spatstat.explore_3.1-0 progressr_0.13.0      
 [46] lmtest_0.9-40          spatstat.random_3.1-4  xfun_0.37             
 [49] globals_0.16.2         timechange_0.2.0       mime_0.12             
 [52] miniUI_0.1.1.1         lifecycle_1.0.3        irlba_2.3.5.1         
 [55] goftest_1.2-3          future_1.32.0          MASS_7.3-58.3         
 [58] zoo_1.8-11             scales_1.2.1           hms_1.1.2             
 [61] promises_1.2.0.1       spatstat.utils_3.0-2   parallel_4.2.1        
 [64] RColorBrewer_1.1-3     yaml_2.3.7             memoise_2.0.1         
 [67] reticulate_1.28        pbapply_1.7-0          gridExtra_2.3         
 [70] stringi_1.7.12         rlang_1.1.0            pkgconfig_2.0.3       
 [73] matrixStats_0.63.0     evaluate_0.20          lattice_0.20-45       
 [76] ROCR_1.0-11            tensor_1.5             labeling_0.4.2        
 [79] htmlwidgets_1.6.2      cowplot_1.1.1          tidyselect_1.2.0      
 [82] parallelly_1.34.0      RcppAnnoy_0.0.20       plyr_1.8.8            
 [85] magrittr_2.0.3         R6_2.5.1               generics_0.1.3        
 [88] DBI_1.1.3              withr_2.5.0            pillar_1.8.1          
 [91] fitdistrplus_1.1-8     survival_3.5-5         abind_1.4-5           
 [94] sp_1.6-0               future.apply_1.10.0    KernSmooth_2.23-20    
 [97] utf8_1.2.3             spatstat.geom_3.1-0    plotly_4.10.1         
[100] tzdb_0.3.0             rmarkdown_2.20         grid_4.2.1            
[103] data.table_1.14.8      digest_0.6.31          xtable_1.8-4          
[106] httpuv_1.6.9           munsell_0.5.0          viridisLite_0.4.1